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ABSTRACT 



Debris disks have been found primarily around intermediate and solar mass stars (spectral types A-K) but rarely around low mass 
M-type stars. We have spatially resolved a debris disk around the remarkable M3-type star GJ 581 hosting multiple planets using deep 
PACS images at 70, 100 and 160 fim as part of the DEBRIS Program on the Herschel Space Observatory. This is the second spatially 
resolved debris disk found around an M-type star, after the one surrounding the young star AU Mic (12 Myr). However, GJ 581 is 
much older (2-8 Gyr), and is X-ray quiet in the ROSAT data. We fit an axisymmetric model of the disk to the three PACS images and 
, found that the best fit model is for a disk extending radially from 25 ± 12 AU to more than 60 AU. Such a cold disk is reminiscent 

QO . of the Kuiper Belt but it surrounds a low mass star (0.3 M Q ) and its fractional dust luminosity Ld us ,/L, of ~ 10~ 4 is much higher. The 

■ inclination limits of the disk found in our analysis make the masses of the planets small enough to ensure the long-term stability of 

the system according to some dynamical simulations. The disk is collisionally dominated down to submicron-sized grains and the 
dust cannot be expelled from the system by radiation or wind pressures because of the low luminosity and low X-ray luminosity of 
GJ 581. We suggest that the correlation between low-mass planets and debris disks recently found for G-type stars also applies to 
M-type stars. Finally, the known planets, of low masses and orbiting within 0.3 AU from the star, cannot dynamically perturb the disk 
over the age of the star, suggesting that an additional planet exists at larger distance that is stirring the disk to replenish the dust. 

Key words, debris disks : circumstellar matter - planetary systems : formation - stars: planetary systems 

X: 

H . 

^.1. Introduction by the evolution of the planetary system in which they form (e.g . 

, „ . „ „ iPetit et alJl200ll iMorhidelli et al.ll2005b iLvkawka et alJl2009h . 

A debris disk around a main sequence star is a collection of small and contain Qbjects whose accretion was stymied b the for . 

bodies left over from the planet formation process. In our Solar mation and mi tion of iant lanets in the tem ^ or sim . 
system, the Asteroid belt and Edgeworth-Kuiper Belt are the two , occmTed tQQ d(Jwl for them tQ , since debrifj 

best known reservoirs of objects that remain from the planet diskfj CQntain a yast number of objects Qn similar ^ 

formation process and range in size from hundreds of kilome- th experience a continua i collisional grinding which produces 

ters in diameter to meter-scale bodies (e.g. | Jewitt et al. | | 2000t and continuall replenishes a population of dust. This dust al- 

| Sheppard&Truiillo || 2006|) . Such reservoirs are highly sculpted 1qws ufj tQ directly detect debrifj diskfj around Qther &ms - m tWQ 

~ ; ~ ; T _ T , ways. The dust is heated by radiation from the central star, and 

Send offprint requests to: J-F. Lestrade, e-mail: . . . ,. . . il , 

jean-francois . lestrade@obspm.fr therefore emits thermal radiation with a te mperature character - 

* Herschel in an ESA space observatory with science instruments ifc of its distance from its host star (e.g | Aumann et al. || l984t 

provided by European-led Principal Investigator consortia and with im- i Greaves et al. || 2005|) . In addition, the smallest grains of dust can 

portant participation by NASA. efficiently scatter the light of the host star (e.g. ISmifh & Terrild 
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Il984t iKalas et al.1 12005). The physical and observation a l prop - 
erties of debris disks were defined by Lagran ge et al.l (12 000). 
and their stud ies were recently reviewed by Wvatt (2008i|) and 
lKrivovl(l2010|) and will event ually place our Solar system in con- 
text dGreaves & Wvattll2010l) . 

Almost a ll debris disks detec t ed by the satell i tes IRAS, ISO 
and S pitzer dBrvden et alJ 120061: ISu et alJ 120061: iTrilling et all 
120081) . HST (Goli mowskietalJl201ll) . and ground-based tele- 
scopes dWvattll2008l) surround A-tvpe and F,G,K-type stars de- 
spite several deep surveys of large samples of M stars conducted 
from mid-IR to submillimeter wavele ngths (Plavchan et alJ 
2005; Lestrade et ail |2006t iGautier et al.1 120071: Lestrade et all 
2009: lAvenhaus et al.ll2012l) . Currently, among the nearby M- 
stars, the only spatially resolved debris disk is around th e very 
young Ml star AU Mic dKalas et alJ l2004t iLiu et alJ |20o4 
iRrist et aD 12005b IWilner et alJl2012[) which has been mod eled 
bv lAugereau & Beusll d2006l) and IStrubbe & Chiang! d2006l) . In 
addition, there are a fe w candidate disks with excesses above 
photospheric level (e .g. ISmith et alJl2006t iLestrade et alll2006t 
iPlavchan et al.l2009l) . Finally, in the cluster NGC2547 (-40 Myr 
old and ~ 433 pc), deep Spitzer MIPS observations have revealed 
1 1 M-stars with 24 pm excesses above photospheric level and 
no excess at 70 pm; these observations have be en interpreted as 
warm dust in debris disks (Forbrich et al. 20081). 

The fact that debris disks are more seldomly observed among 
M-stars than around higher-mass stars seems surprising at first, 
since all spectral types have similar detection rates of proto- 
planetary disks in the earlier stage of their evolution, accord- 
ing to observati ons of low density clusters like Taurus-Auriga 
and p Oph (e.g.. lAndrews & Williams! 120051) . However, in high 
density clusters like Orion, external photoevaporation by intense 
FUV radiation field can severely limit the production of plan- 
etesimals around low mass M -stars on a timescale shorter than 
~10 Myr dAdams et al.l l2004). Another hazard for M-stars dur- 
ing the first ~100 Myr is close stellar flybys, when co-eval stars 
are still in the expanding cluster of their birth and strongly in- 
teracting with each other. During these early close stellar flybys, 
planetesimals are stripped from disks, and this is more severe for 
disks around low mass stars in high stellar density c lusters like 
Orion accordi ng to simula t ions ([L estrade et al. 201 1). 

Recently, Wva tt et all d2012l) have found evidence of the 
prevalence of debris disks in low-mass planetary systems (also 
Moro-Martin et al. in prep) and suggest that this correlation 
could arise because such planetary systems are dynamically sta- 
ble over Gyr timescales. Recent observations show that low- 
mass planet s are more abundan t among M-stars than around the 
other stars (Bonfils et al .1 1201 It lHowardetaLll2012l) . Hence, if 
the correlation between debris disks and low-mass planets for 
G-stars applies to M-stars, then debris disks should be relatively 
common around them, in contrast to a paucity of detections. 

However, debris disks around M-stars are harder to detect 
than around more massive stars at the same distance simply be- 
cause they are less luminous, meaning that the dust within expe- 
riences significantly less heating. Therefore, to detect the same 
disk around a later type star requires deeper observations. M-star 
debris disks may also be less detectable because additional grain 
removal processes are operating. For example, a physical pecu- 
larity of M-stars is that they are structurally different from solar- 
type stars. Their interiors have deep convective zones - fully 
convective for M3 spectral type and later - that produce strong 
coronal magnetic fiel ds responsible for th eir optical/radio flares 
and X-ray emission dHawlev et alJl2000h . T his activity gener- 
ates a lso stellar winds of energetic particles ( Wargeli n & Drake] 
120011) which might dominate the circumstellar grain removal 



proce sses for a large fraction of the star lifetime dPlavchan et al.1 
120051) . 

This paper describes observations carried out as part of 
the Key Program DEBRIS (Disc Emission via a Bias-free 
Reconnaissanc e in the Infrared/S ub-mm) on the Herschel Space 
Observatory dPilbratt et al.ll2010l) . DEBRIS is an unbiased flux- 
limited survey to search for dust emission at A = 100 and 
160//m toward the nearest ~89 stars of eac h spectral type 
A,F,G ,K,M as evidence of debris disks (see Mat thews et al.l 
(120 1 Oh and iPhillips et all d2010l) for the sample description). 
For selected targets, complementary Herschel observations at 
70, 250, 350, 500 pm were also conducted. The first results 
of this program have already shown that these observations 
can detect disks down to much fainter levels than previously 
achieved, and mor eover can spatially resolve debris disks at 
far-IR wavelengths dMatthews et"aj]l2010l:IChurcher et afli 20 1 1 ; 
Kennedy et"aT]|2012at IWvatt et alj|2012t Kennedy et alJl2012bT 
iBooth et aUl2012l : iBroekhoven-Fiene et al.ll2012l) . 

As part of this survey, we have spatially resolved a disk 
around the M3 spectral type star GJ 581 at A=70, 100, and 
160 pm. Hence, this is the second resolved debris disk around 
an M-star, but, in contrast to the sta r AU Mic which is young 
(12 Myr, IZuckerman & Songl 12001) . GJ 581 is old (2-8 Gyr, 
see § [3j. Also, GJ 581 is surrounded by at least four low mass 
planets with minimum masses of 1.9, 15.6, 5.4, and 7.1 M e , 
orbital radii of 0.03, 0.04, 0.07, and 0.22 AU, and eccentrici- 
ties between 0.0 and 032, detected by radia l velocity me asure- 
ments dBonfilset alJ l2005t lUdrv et alj|2007l: iMavor et alj[2009t 
iForveille et alj|201lb . All these planets are within the tidal lock 
region of this M3 spectral type star (<0.25 AU) and hence are 
expected to be synchronousl y rotating and potentially undergo- 
ing at mospheric instabilities dWordsworth et alj|201 UlKite et al.l 
1201 ll) . Planets GJ 581c a nd d are near and in the conventionally 
defined Habitable Zone dSelsis et al.l [2007), respectively. The 
presen ce of one or two additional planets in the system is d e- 
bated dVogt et al.ll20Tol: IForveille e"tai1l201 UlVogt et al.ll2012l) . 

In this paper, we describe the Herschel observations of 
GJ 581 as well as archival MIPS and IRS data from Spitzer, and 
NICMOS data from HST in § [2] The stellar parameters of GJ 5 8 1 
used are in § [3] Reconnaissance of a cold debris disk around 
G J581 in the three PACS images at 70, 100, and 160 pm and 
in the presence of background sources contaminating the field 
is described in § [4] Modeling of these images to determine the 
spatial distribution of the emitting dust is described in § [5] The 
spectral energy distribution SED including the IRS spectrum of 
GJ 581 and modeling of a hypothetical second component of 
warm dust are described in § [6] An upper limit on the bright- 
ness of scattered light using the NICMOS image is discussed in 
§ |7] Physical conditions in the disk and its relationship with the 
planetary system around GJ 581 are discussed in §[8] 



2. Observations 

2.1. Herschel 

GJ 581 was initially observed w ith PACS (Photodet ector and 
Array Camera & Spectrometer, Pogli tsch et al.l (1201 Oh ) on 1 1 
August 2010 using the standard DEBRIS observing strategy, and 
a resolved disk was tentatively detected at 100 and 160 pm. We 
then acquired deeper PACS images at 100 and 160 pm on 29 July 
2011, a PACS image at 70jum (and 160yL/m) on 1 August 2011, 
as well as SPIRE images at 250, 350 and 500 pm on 30 January 
201 1. These observations are summarised in Table [TJ 



2 



J.-F. Lestrade et al.: Resolved Disk Around GJ 581 



Table 1. Herschel observations of GJ 581. 



Obsld 


Date 


Instrument 


Integration 


1342202568 


11 August 2010 


PACS 100/160 


890s 


1342213474 


30 January 2011 


SPIRE 250/350/500 


185s 


1342224948 


29 July 2011 


PACS 100/160 


7190s 


1342225104 


1 August 2011 


PACS 70/160 


3936s 



2.1.1. PACS 

The PACS observations used the mini-scan map mode with eight 
legs of a 3' length, with a 4" separation between legs in a single 
scan direction at a rate of 20 arcsec s , and two scan directions 
(70° and 110°). These data were reduced using the Herschel 
Interactive Processing Environment HIPE dQt3[2010) version 7 
and implement version FM6 of the flux calibration. The data 
were pre-filtered to remove low-frequency (1//) noise using a 
box-car filter with a width of 66 arcsec at 70 and 100 //m and 
102 arcsec at 160 /mi. This data filtering results in the source 
flux densit y being underestimated by ~ 20 + 5% as discussed 
in detail bv lKennedv et al.1 d2012al) . Maps were made from these 
filtered timelines using the photProject task in HIPE. 

The pixel scales in the images presented in Fig Q] were set to 
1 arcsec at 70 and 100 yum, and 2 arcsec at 160 fim, i.e., smaller 
than the natural pixel scales. This enhanced sampling is possible 
because of the high level of redundancy provided by the scan 
map mode used but it comes at the cost of correlated noise be- 
tween neighbouring pixels. We have also made images with the 
natural pixel scales of 3.2 arcsec at 70 and 100 /im, and 6.4 arc- 
sec at 160 yum to evaluate the impact on the parameter estimation 
in our modeling. The noise rms for the images with the natu- 
ral pixel are 0.47mJy/5.6"beam at 70 /-im, 0.48mJy/6.7"beam at 
100 jum, and 0.77mJy/l 1.4"beam at 160 //m. 



2.1.2. SPIRE 

Follow-up observations were taken on 30 Januar y 2011 with 
SPIRE (Spectral & Photometric Imaging REceiver, Griffin ~et~ai] 
(2010)) using the small-map mode, resulting in simultaneous 
250, 350 and 500 //m images. The data were reduced using 
HIPE (version 7.0 build 1931), adopting the natural pixel scale 
of 6, 10, 14 arcsec at 250, 350 and 500 yum respectively. The 
noise rms are 6.1 mJy/18.2"beam, 7.9 mJy/24.9"beam, and 
8.3 mJy/36.3"beam at 250, 350 and 500 /urn, respectively, and 
the image at 250 yum is shown in § 14.41 



2.2. Ancillary data 
2.2.1. Spitzer 

MIPS 70 jum observations of GJ 581 (AOR 22317568) were 
taken on 21 August 2007 (no 24 /jm MIPS were taken) and 
a small measured ex cess, with the signifi cance xio = (F°q S ~ 
F* Q )/o-7o = 3.6 in iKospal et al.1 d2009l) and xio = 2.2 in 
Brvd en et al.l (120091) with the same data, forms a tentative 
discovery. We re-reduced the archival data usin g an updated 
pipelin e and the flux calibration summarized in Gordo n et al.l 
(2007) providing the new flux density 20.0+5.3 mJy by PSF 
fitting to the image at the effective wavelength of 71.42 yum 
(color correction for Trf„ Jf =40 K applied : 0.89 Q). The uncer- 



tainty includes both statistical and calibration uncertainties and 
the corresponding excess ratio xio is 2.7 with our estimate of 
the photospheric flux density of 5.6 mJy at 71.42 fim (see §[6}. 
Photospheric flux densities predicted for late type stars (K and 
M) by the Kurucz or Next Gen models have bee n shown to be 
overestimated in the mid- IR by as much as 5-10% dGautier et al.l 
12007b lLawler et al.l l2009). Hence, this excess can be treated as a 
lower limit. 

The IRS observations of GJ 581 (AOR22290432) were taken 
on 3 1 August 2007, an d th e details of the data reduction are in 
iBeichman et al.l (120061) and lDodson-Robinson et al. 

2.2.2. HST/NICMOS 

GJ 581 was directly imaged with HST/NICMOS on 6 May 
1998 (GO-7894; PI Todd Henry). The NIC MOS data and the 
ove rall observing program a re described in iRrist et aTl (jl998) 
and iGolimowski et al.r(l2004l) . We reanalyzed the F110W data 
for GJ 581 consisting of 128 seconds of cumulative inte- 
gration on the NIC2 camera (0.076"/pixel, 256x256 pixels). 
Target stars were not placed behind the occulting spot, near- 
contemporaneous observations of PSF reference stars were not 
made, and multiple telescope roll angles were not employed. 
Therefore the observations were not optimized for high-contrast 
imaging of low surface brightness circumstellar nebulosity. 
Nevertheless, we subtract the GJ 581 point-spread-function 
(PSF) using observations of LHS 1876 (GJ 250B) made on 24 
March 1998 as part of the same scientific program, and in so 
doing, set constraints on the scattered light disk brightness as 
discussed in ^8] PSF subtraction techniques, including a discus- 
sion of scattered light artifacts and other spuriou s features, are 
described in greater 



ht artiracts and other spurious 
detailed bv lKrist et all (Il998l) . 



1 http://irsa.ipac.caltech.edu/data/SPITZER/docs/mips/ 
mipsinstrumenthandbook/5 1/ 



3. Stellar parameters of GJ 581 

GJ 581 (HIP 74995) lies relatively nearby (6.338 + 0.071 pc ; 
Philli ps et ail (1201 0l)) an d is classified as a star of spectral type 
M3.0 dReid et al.lll995h . Recent CHARA interferometric mea- 
surements of its physical radius (0.299 ± 0.01QR o ) imply an ef- 
fective surface temperature of T e ff = 3498 ± 56K, a bolometric 
luminosity of 0.0 1205 + O.OOO24L and a stellar mass of O.28M 
(Ivon Braun et al.ll201 ll) . 

A variety of different techniques have been discussed in the 
literature as a means to determine the age of GJ 581, includ- 
ing kinematics, magnetic activity (X-ray observations), chromo- 
spheri c activity, stellar color, metallicity and rotation. Leggett 
(1992) finds that the galactic velocities of GJ 581 are interme- 
diate bet ween those typical of the young and old galactic disk 
M-stars. [Bonfils et al.l (12005) conclude that the low limit on its 
X-ray emission, the low v sini and the weak Call H and K emis- 
sion, taken altoget her, suggest that GJ 581 is at least 2 Gyr old. 
ISelsis et al.1 (2007) established an L JLt, i versus age relation for 
M- K- G- spectral type stars to es timate that the age of G J 581 
could be around 7 Gyr. Recently, lEngle & Guinanl (1201 ll) have 
established an age-rotation period relation for M-stars and deter- 
mined an age of 5.7 + 0.8 Gyr for GJ 581. Clearly, GJ 581 is an 
old star well above 1 Gyr. 

High contrast imaging for GJ 58 1 has revealed n o companion 
of ~7 Jupiter masses or higher between 3-30 AU dTanner et al.l 
2010). In addition, the limits provided by the HARPS radial ve- 
locity measurements exclude planets that are more massive than 
Jupite r with semimajor axes inside 6 AU (Fig 13 in IBonfils et al.l 
120111) . 
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All the parameters of GJ 581 are summarized in Table [2] 



Table 2. Stellar parameters of GJ 581. 



Parameters 


values 


references 


R.-A. ICRS(2000) 


15 H 19 m 27 509 s 


Hete et al (20001 


Dec. ICRS(2000) 


-07°43' 19.44" 


fx a cosS 


-1.228"/yr 






-0.098"/yr 


» 


Galactic longitude 


354.08° 




Galactic latitude 


+40.01° 


» 


Distance 


6.338 ±0.071 pc 


Phillips et al. (20101 


Spectral type 


M3.0 


Reid et al. (1995) 


Radius 


0.299 ± 0.010 Rq 


von Braun et al. (201 1) 


Mass 


0.28 M 




Bolometric lum. 


1.22 ±0.02 10~ 2 L o 




Effective temp. 


3498 ± 56 K 




Metallicity [Fe/H] 


-0.25 


Bonfils et al. (20051 


v sini 


< 2.1 km/s 


Delfosse et al. (19981 


Rotation period 


93.2 ± 1 days 


Vogtetal. (20101 


LogLx (ergs/s) 


< 26.44 


Schmitt et al. (19951 


Age 


2 - 8 Gyr 


see §[3] 



4. Herschel images of GJ 581 

In FigQ] we present our deep PACS images of GJ 581 cropped 
to the region ±50" from the star. At each wavelength, the main 
emission is close to the star position (image center) and the sur- 
rounding field is contaminated by several other sources, detected 
barely at 70 fim, significantly at 100 fim, and prominently at 
160 fim, as is expected for submillimeter galaxies in the back- 
ground. The central emission is suggestive of a spatially resolved 
debris disk not fully separated from the background source to 
the N W. In the next subsections, we analyse quantitatively these 
PACS images to verify this view. 

4.1. Radial profiles of the emission in the PACS images 

First, in Fig[2j we present the radial profiles of the emission at the 
three wavelengths by computing the mean brightness in ellipical 
annuli centered on the peak of the main emission. At the three 
wavelengths, these radial profiles show that the emission is ex- 
tended about this peak when they are compared to the Gaussian 
profiles of a hypothetical point source in the background. 

We computed these profiles after subtracting from each im- 
age a PSF scaled to the photosphere flux density (S photosphere = 
5.8,2.8, 1.1 mJy at 70, 100 and 160 fim, respectively, see ©. 
The emission peak position was found to be less than 2 arcsec 
from the star in each image, consistent with the pointing accu- 
racy of the Herschel telescope 0. We have tested elliptical annuli 
at PA =120° and with inclinations 0° (circular), 25°, 40° and 
75°, anticipating that the disk may not be in the plane of the 
sky. We found that all these radial profiles were more extended 
than the Gaussians at the three wavelengths, and slightly more 
for the inclination of 40°. The imprint of the N W source at the 
radial distance of 11 arcsec can be seen in these profiles at 100 
and 160 /im. We have tested the method by computing the ra- 
dial profile of the emission of the South East background source 
in the same way. We satisfactorily found that its profile matches 
within lcr the Gaussian expected for a point source. 

2 http://herschel.esac.esa.int/twiki/bin/view/Public/SummaryPointing 



The extended emission revealedby these profiles can also 
be seen directly in the photosphere-subtracted images of Fig Q] 
(middle column), most prominently at 70 fim because the pho- 
tospheric flux density is highest at this wavelength. 

4.2. Gaussian source fits 

As a second approach to verify that the central emission is more 
extended than the PSF, we fit an elliptical 2D Gaussian to each 
photosphere-subtracted image with masking applied to the po- 
sition of the N W source at 100 and 160 fim. At 70 /im, we 
found FWHM of the minor and major axes of 9.5 ± 0.7" and 

10.7 + 1.1"; a position angle ofl20 ± 10°; and a flux density of 

18.8 ± 1.4 mJy after adding back in the photospheric contribu- 
tion. At 100 //m, we found, respectively, values of 9.9 + 1.0", 
13.3 ± 0.5", 120 ± 10°, and 21 .1 ± 1 .5 mJy for the minor and ma- 
jor axes FWHM, the position angle, and the 100 fim flux density, 
after masking the N W source (all pixels in a 12" x 12" square 
centered at -9" and +6" from the star). Given that the FWHM 
of the PACS PSFs are 5.6" and 6.8" at these wavelengths, we 
conclude that the emission is significantly extended, as already 
shown by the radial profiles, and is elongated at PA ~ 120°. 

The ICRS coordinates of the 70 fim Gaussian peak in this fit 
are 15 h 10 m 25.905+0.05 8 and -7° 43' 22.79 ± 0.7" and differ 
only by 0.3" from the adjusted position of the 100 /im Gaussian 
peak. These coordinates differ by +0.66" and - 1 .48" from the 
right ascension and declination of the star GJ 581 predicted with 
the Hipparcos astrometric parameters (Table [2J. These differ- 
ences are consistent with the lcr pointing accuracy of 2" for the 
Herschel telescope. We conclude that the main emission in the 
PACS images at 70 and 100 /mi is centered on the star position 
within pointing uncertainty. 

The 160 fim image is the product of the coaddition of two 
images taken independently with the PACS 70/160 and PACS 
100/160 instruments only a few days apart in 2011 (Table [TJ. 
The registration of these two images were facilitated by the neg- 
ligible displacement due to proper motion over the short lapse of 
time and by the fortuitously small difference of 0.3" between the 
pointing positions of the two instruments as found above. Our 
first 160 fim image in 2010 was not coadded because no fea- 
ture in the image could be used to check the registration since its 
signal-to-noise ratio was V12.5 times lower. Although this first 
image was crucial in our decision to observe deeper, its use or 
not is inconsequential for our analysis. For the fit at 160 fim, we 
fixed the position of the Gaussian to the coordinates determined 
at 70 fim. This was necessary because the large mask (16" x 16" 
square) used for the N W source affected the independent deter- 
mination of this position. The best fit 2D Gaussian parameters 
were minor and major axes of 12.8 + 1.5" and 21.5 ± 2"; a posi- 
tion angle of 125 + 10°; and a flux density of 22.1 ±5.0 mJy, after 
adding back in the photospheric contribution. Given the PACS 
PSF at 160 fim of 11.4", this indicates that the emission is ex- 
tended at this wavelength as well. 

The disk inclinations resulting from the ratios of the minor 
and major axes determined above, and corrected quadratically 
for the convolved PSFs, are : 33 ± 17° (0° is face-on), 54±6°, and 
71 ± 7° at 70, 100 and 160 fim respectively. Although scattered, 
these values are statistically consistent, since they are within 
1.5cr from their weighted mean (59°), and indicate an inclined 
disk which has implications for the masses of the planets of the 
system as discussed in § 18.31 

We note that the three major axes above, corrected quadrati- 
cally for the convolved PSF, are very closely proportional to the 
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Fig. 1. PACS images of GJ 581 cropped to +50 arcsec from the star, at 70, 100 and 160 //m from top to bottom, respectively. In the 
left-hand column, the raw images show that the main emission is centrally located about the star position (image center) and that 
there are several point-sources in the field, barely detected at 70 yum, significantly at 100 yt/m, and more prominently at 160 pm, as 
expected for submillimeter background galaxies. In §j4] we show that the main emission is extended and centered on the star position, 
as expected for a debris disk, and mingles with a background source ~1 1 arcsec toward the North West. The panels in the middle 
column are the photosphere-subtracted images. The panels in the right-hand column show the best subtraction (lowest residuals) 
of a two-point source model, which assumes that there is no debris disk around GJ 581 but an extra background source located 
exactly behind the star in addition to the N W source. This model is rejected because of the systematic residuals left, indicative of 
an extended structure, especially at 70 and 100 pm. At each wavelength, the contours levels of the three images are the same and 
correspond to l,2,3,9,15o- (<x = 0.0135 mJy/l"pixel) at 70 /mx, l,2,3,6,9,12,15cr (cr = 0.0094 mJy/l"pixel) at 100 pm, and 

I, 2,3,5,7,9,1 lo"o (o"o = 0.0251 mJy/2"pixel) at 160 pm. The coordinates of the image center provided in the labels correspond to 
the star position at epoch of observation (July 29th - August 1st 2011). The hatched circles are the beam FWHMs : 5.6, 6.8, and 

I I. 4 arcsec at 70, 100, and 160 /mi, respectively. 



wavelength, suggesting that the disk is radially broad since emis- 
sion at longer wavelength probes colder dust, more distant from 
the central star. 



The flux densities from our fits above have been scaled up to 
account for the flux removed by the data filtering during the re- 
construction of the images; the correction factors are 16, 19 and 
21% with an uncertainty of 5% estimated for point sources in 
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Fig. 2. Radial profiles of the mean brightness of the photosphere-subtracted images. Each point of these curves was computed as the 
mean of the emission in an elliptical annulus centered on the peak of the main emission close to the image center. They correspond 
to an axisymmetric disk model inclined to 40° relative to the plane of the sky (see text). They are one pixel wide, i.e. 1 arcsec in the 
70 and 100 /mi images and 2 arcsec in the 160 /mi image. The Gaussian with the FWHM of the beam is the profile expected for a 
hypothetical point source in the background aligned by chance with the star and shown for comparison. 



the DEBRIS survey bv lKennedv et alj d2012al) . The uncertainty 
of each PACS flux density determined above is based on the 
quadratic sum of the statistical uncertainty in our Gaussian fit, 
the absolute flux calibration accuracy of 3% at 70 and 100 /mi, 
and 5% at 160 /an, provided by the Herschel project^, and the 
5% uncertainty of the correction factor for the data filtering. The 
uncertainty of the flux density at 160 fim is formally 2 mJy with 
this calculation but our fit depends to some degree on the posi- 
tion of the mask applied for the N W background source. So we 
have increased this uncertainty to 5 mJy at this wavelength based 
on several fits with different masks. 

4.3. Considering the superposition of two point-sources 

As a final test, we consider the possibility that the extended emis- 
sion in the central part of the image could be caused by the su- 
perposition of two backgound sources instead of a disk. To this 
end, we subtracted two PSFs from each PACS image (in addi- 
tion to subtracting the phostospheric emission), adjusting their 
flux densities in order to remove as much emission as possible. 
The first PSF was located at the position of the N W source, i.e. 
at -9" and +6" from the star, and the second was tested at six 
locations ; the star position itself, as well as a half FWHM to the 
North, South, East and West of the star, and half way between 
the star and the N W source. The lowest residuals were found 
after removing 1.4, 2.9 and 6.6 mJy for the first PSF at 70, 100 
and 160 /mi, respectively, and 4.9, 6.9 and 9.4 mJy for the second 
PSF at the star position at 70, 100 and 160 /mi, respectively. Note 
that these latter flux densities are free of the photosphere contri- 
butions estimated in § 16.11 Despite this removal process, there 
is still significant structure left in the residual images shown as 
the two-point source subtracted images in Fig [T] (right-hand col- 
umn). This structure can be best explained as resulting from the 
extended emission of the disk incompletely removed by this pro- 
cess. Hence, we conclude that this test rejects the possibility that 
the superposition of two background point sources can be re- 
sponsible for the central emission. 

We elaborate further by discussing the probability to find 
such contaminant sources in the field and their spectra. First, 
the probability to have one background source stronger than 
6.6 mJy at 160 //m within 11" from the star is 18 %, and to 



3 http://herschel.esac.esa.int/twiki/pub/Public/PacsCalibrationWeb/ 



have two is only 1.8% by using the Poisson probability distri- 
bution with the mean source surface density N(S > 6.6mJy) ~ 
4000 sources/deg 2 a t 160^/m pr ovided by the Hersche l PEP sur- 
vey (iBerta et al.l2011l) (see also lSibthorpe et al.l(l2012l) ). Second, 
the spectra of these test sources removed from the images may 
or may not be physical. The flux densities removed at the 
N W source position (1.4, 2.9 and 6.6 mJy at 70, 100 and 160 /mi 
respectively) are consistent with the spectrum of a galaxy at 
z > 1 .5 according to Fig 4 o f lBlainetal] (12002) valid for a typi- 
cal high-z galaxy enshrouded in dust (~ 38 K, L ~ 5 x 10 12 L o ) 
and radiating in the far-IR and submm. However, the flux densi- 
ties removed at the star position (4.9, 6.9 and 9.4 mJy at 70, 100 
and 160 //m respectively) above the photospheric levels make the 
ratio S ioo/jm/S70/jm lower than expected for a galactic spectrum 
according to that work. 

4.4. SPIRE images 

The SPIRE image at 250 /mi in Fig [3] shows an elongated struc- 
ture at PA ~ 120° which has two peaks at the 2cr level that match 
the positions of the star and the N W source. This structure is 
also extended to the S E. Although this structure is of low sta- 
tistical significance, it is consistent with the emission detected at 
the PACS wavelengths. The two other SPIRE images, at 350 and 
500 /mi, are dominated by noise and no reliable structure can 
be recognized. A 2D Gaussian could not satisfactorily model the 
emission at 250 /mi, and so we carried out photometry with an 
aperture of 36 arcsec and measured a flux density of 24 + 6 mJy. 
This flux density is considered an upper limit for the disk since 
its emission is blended with that of the N W source. 

4.5. Spitzer MIPS image at 70 pm 

A Spitzer MIPS image at 70 /mi was made on 21 August 2007, 
four years before our PACS image (1 August 2011). We fit a 
Gaussian with the FWHM of 19 arcsec (70 /on MIPS beam) 
to the emission of the MIPS image and found the coordinates 
of the Gaussian peak to be at a = 15 h 10 m 26.08±0.17 s and 
5 = -7° 43' 23.2 + 2.5" in the ICRS system. The relatively 
high noise of the MIPS image did not permit a solution for 
the FWHM. The differences in the coordinates compared with 
our 70 /mi PACS position given in § 14.21 are Aa = -2.56" 
and A6 = +0.45" (PACS minus MIPS) with an uncertainty of 
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Fig. 3. SPIRE image of GJ581 at 250 pm. Pixel size is 6" and 
the contour levels are 1 and 2cr with cr- 6.1 mJy/18.2"beam. The 
star symbol is the position of GJ 581 at the date of observation. 

3.3" when combining the astrometric uncertainty of the MIPS 
Gaussian fit (2.5" in both coordinates) and the pointing accura- 
cies of 2" for Herschel and of 1" for SpitzerQ The predicted 
displacement of GJ 581 between the two epochs of observations 
is Aa = -4.98" and AS = -0.41" computed with the proper mo- 
tion in Table [2] Hence, the coordinate differences of the 70 pm 
emission measured between the two epochs are compatible with 
this prediction but the star has not moved sufficiently between 
these epochs for us to confirm that the 70micron emission is co- 
moving with the star at a statistically significant level. 

5. Modeling the PACS images of GJ 581 

5.1. Parametrized model 

We fit a parametrized model of the disk to the PACS images. 
The model is axisymmetric and truncated by the inner radius r m 
and the outer radius r out which are free parameters in our fit. Its 
dust emission is optically thin, and the flux density from each 
element (k, I) of the grid covering the disk is 

AS kJ = eB v (T(r kJ )) . . Aa/d 2 , (1) 

where B v {T{rkj)) is the Planck function that depends on the grain 
temperature T(fkj) at the radial distance from the star, Aa 
is the area of the element in the grid, d is the dist ance to the 
star, and T. p is the coefficient of the power-law (e.g. Wva tt et all 
1999) B To fit individuaUy each PACS image in § E2 we set 
the factor e to unity in eq (1). To fit simultaneously the three 
PACS images at A = 70 pm, 100 pm, and 160 pm in § 15.31 we 
implement a grey body effect (e.g. lDent et al.ll2000l) by setting e 
to unity if A < Aq, and e = 1 .0 x (Ao/A~f if A > Ao, where Aq and 
B(> 0) are free parameters in our fit. 

In our model, no assumption is made for the size distribution 
of the grains, their mineralogical composition and porosity. The 
thermal structure of the disk is taken as 

T(r) = f T . T BB (r), (2) 

where fx is a scaling factor applied to the black body tempera- 
ture T BB (r) = 278 . (L/L ) - 25 . (r/lAUy 05 (K), and is a free pa- 

4 http://irsa.ipac.caltech.edu/data/SPITZER/docs/spitzermission/ 
missionoverview/ spitzertelescopehandbook/ 1 2/ 

5 The re is a different bu t equivalent derivation of the flux density 
given in IZuckermarJd200lh based on the surface emittance nB v . 



rameter in our fit. Here we illustrate how to interpret this param- 
eter using a simplified model of the absorption and reemission 
of the starlight by the grains. For dust with the absorption and 
emission efficiencies Q a i, s and Q, a straightforward derivation 
shows that fx = (Qabs/Q) 1 ^ 4 for a grain at thermal equilibrium, 
ignoring that these efficiencies should be averaged over the spec- 
trum of incoming and outgoing radiation, and integrated over the 
dust size distribution. If we assume that grains larger than 1 pm 
absorb starlight with the efficiency gafes - 1 an d reemit it at 
longer wavelengths at the lower efficiency Q, then the simple in- 
terpretation is that fx is directly related to the emission efficiency 
through the relationship fx = Q . 

The term "Lpf* in eq (1) is the emitting cross-sectional area of 
the grains per unit area of the disk surface. These grains are spa- 
tially distributed according to a radial profile taken as the power- 
law "Lpr®, and their total cross-sectional area is A. Since these 
grains reemit with the efficiency Q, their total emission is pro- 
portional to Q . A = JT 2nrdrL p r a . Hence, if a + -2, the 
coefficient of the power-law is 

Z p = f T 4 . A . (a + 2) / (2n (C? - r^ 2 )), (3) 

and the total cross-sectional area A and the power-law index a 
are free parameters in our fit 0. Grains smaller than 1 pm are 
not important becau se they emit so ineffici ently that their flux 
density is neglig ible (lBonsor&Wvattll2010h . 

The total cross-sectional area A can be converted into mass 
assuming a size distribution and a mass density p for the ma- 
terial. Adopting the standard size distribution n(D) oc D~ 3 5 for 
spherical particles of diameter D between D m ,„ and D max , the 
corresponding mass is 

m d = (2/3) .A. p. y/D mi „ . y/D max . (4) 

This model is complemented with a point source photo- 
sphere centered on the image by two free parameters (coordi- 
nates x c and y c ) and having the flux densities estimated from the 
Next Gen stellar atmosphere model in §|6]but lowered by ~ 20% 
because of the data filtering used to reconstruct the images as 
already mentioned in § 12.1.11 This model is projected onto the 
sky with the inclination i and the node orientation CI, and finally 
convolved by the telescope PSF provided by the PACS images of 
the reference star a Boo at 70 pm, 100 pm and 160 pm. Hence, 
our model has 9 free parameters (r,„, r out , a, fx, A, i, £2, x c ,y c ). 

The best fit is found by minimizing 

2 V V / ~ Ck,l\ 2 ,„ 
k I 

computed with the residuals between the image (O) and the 
model (C) over all the pixels of the image, and assuming 
the same measurement uncertainty <tq for all the pixels. To 
compute the model, we used a grid on the sky which has a 
resolution of 0.5" at 100 pm and 70 ^m, and 1" at 160 pm, 
i.e. twice as fine as the pixel size of the images of Fig [TJ and 
has dimensions 128 x 128 at 70 and 100 pm, and 64 x 64 
at 160 pm. These dimensions can accommodate the largest 
disk model tested (r OI „=150AU) extended by twice the beam 
FWHM. This sky grid is the same for all the models tested so 
that the number of degrees of freedom v is the same for all 
of them. To use the xl probability distribution to discriminate 
between them, we carefully estimate the a priori uncertainty cr () 

6 if a = -2, E p = f- 4 .A I {2n (ln(r m „) - ln(r m ))) 
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Table 3. Best fit models of the disk for each individual image and for the combined fit. 



image(s) 


X y 
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A " 
(AU 2 ) 


' inner 

(AU) 


r c 

(AU) 
lower (formal) 


a 




ft 


i d 

C) 


O 




70 pm 


1.04 


408 


0.8 ± 0.5 


25+20 

-14 


> 40 


(62) 


— 2 < a 


< 


4 +0 - 5 


< 60 


120 ± 


20 


100 //m 


0.99 


408 


2.2 ± 1.2 


31 +22 

■"-14 


> 55 


(100) 


-2 < a 


<0 


3.5!f 


< 60 


110 ± 


20 


160 //m 


0.97 


90 


1.5 ±0.7 




> 90 


(145) 


-2<a 


<0 


4.5!f 


40 < i < 80 


120 ± 


20 


Combined 


1.03 


922 


2.3 ± 1.1 


25!|? 


> 60 


(110) 


-2<a 


< 


3.5!?- 5 


30 < / < 70 


120 ± 


20 



" : based on the number of natural pixels (the natural pixel is ~ 3 times larger than the pixel in the images shown). 
* : A is the total cross-sectional area of all the grains (see § 15. It 

c : lower limit (see text) and, in parenthesis, the formal value corresponding to the minimum^ given. 
d : i = 0° is face on. 
e : Q. > is E of N 



by computing the noise rms over the image limited to the sky 
grid dimensions and after excluding the central part (r < 25") 
where the disk emission is located. For GJ 581, the noise rms, 
cr , is 0.0 135mJy/l "pixel at 70 pm, 0.0094mJy/l "pixel 
at 100 /mi, and 0.025 lmJy/2"pixel at 160 pm, corre- 
sponding to 0.49mJy/5.6"beam, 0.50mJy/6.7"beam, and 
0.81mJy/l 1.4"beam, respectively, by using the beam area 
7rxFWHM 2 /41n2. 

Finally, we used the SPIRE 3<x upper limits as a constraint to 
reject any model with a flux density of the dust emission larger 
than 24 mJy at 250 /mi. 

5.2. Fits of individual PACS images 

First, we searched for the best fit model for each individual im- 
age. The ranges of the model parameters tested were: A of 1 to 
20 AU 2 ; r m of 3. 1 to 80.0 AU; r ou , of r, to 150.0 AU; a of -3.0 to 
0.0; f T of 1.0 to 6.0; of 0.0 to 90.0° (0.0° is a disk seen face-on); 
and Q of 0.0 to 180.0° (0.0° is North and increasing Q is East). 
The range for the index a was chosen to cover possibilities such 
as grains being blown out of the system by radiation pressure 
(a = — 1), and having a distribution akin to that of the Minimum 
Mass Solar Nebula (MMSN) (a = -1.5), as d iscussed in the 
model ing of the disk around the A-star /3 Leo bv lChurcher et alJ 
d2011l) . The two background sources, N W and S E of GJ 581, 
were masked as shown in the residual maps of Fig [4] Roughly 
50 million models were tested in our search for the best fit. 

We found the reduced^ 2 = 1.04, 0.99, and 0.97 for the best 
fits to the three images at 70, 100, and 160 /mi, respectively. The 
numbers of degrees of freedom are v = 408 at 70 and 100 /mi, 
and v = 90 at 160 pm. These reduced^ 2 indicate noise-like post- 
fit residuals according to^ 2 -statistics. The residual maps in Fig|4] 
do not show any systematic residuals as expected in these con- 
ditions. We stress that the uncertainty <x used for eq (5) is the 
value determined a priori and is not purposefully tweaked a pos- 
teriori to make the reduced^ 2 close to unity. 

The best fit values of the parameters are in Table [3] There 
are significant correlations between parameters, especially be- 
tween A, a, r outer and fj, as we found by inspecting the two- 
dimensional projections of the x\ hypersurface, e.g. Fig [5] for 
the pair A and fj. To estimate the parameter uncertainties in 
these conditions, we have determined the lower and upper limits 
around the best fit value of each parameter that correspond to the 



fits in which the reduced^ 2 are increased to 1.12 and 1.25 with 
all the other parameters freely adjusted. These thresholds cor- 
respond to a probability of 5% in ^-statistics that the reduced 
xl of pure noise exceeds 1.12 and 1.25 for the number of de- 
grees of freedom v = 408 and 90, respectively. This is a stan- 
dard criterium in fitting procedures. We have also inspected the 
corresponding residual maps and noticed nascent systematics for 
these degraded fits as expected. For the outer radius r out , only the 
lower limits and the best fit values of the fits are given in Table 
[3] because the upper limits are not well constrained since any 
distant dust becomes very cold, even accounting for the incident 
inters tellar radiation field as a source of heating (Lestrade et al. 
2009, Fig Al). The resulting range for r,„ and r OM does not per- 
mit a conclusive estimate of the radial breadth of the GJ 581 
disk. 



5.3. Combined fit of the three PACS images 

In order to consolidate these results and to break correlations 
between parameters, we combined the three PACS images at A = 
70, 100, and 160 pm in a single fit by setting the factor e of 
eq. (1) to unity if A < A and to e = 1.0 X (Ao/Af if A > A , 
in order to implement a grey body effect. Our search covered 
successively the combinations of (3 — 0.0, 0.5, 1.0, 1.5, 2.0 and 
Ao = 70 pm, 85 /mi, 130 /mi. The ranges of the other model 
parameters tested were the same as for the individual images in 
§ 15.21 The two background sources were masked. Note that the 
160 pm image with 4 times fewer pixels has a lower weight than 
the two other images in this combined fit. 

The best fit has the reduced^ 2 of 1.03 (v= 922) for /3 = 0, 
indicating formally no grey body effect for Aq < 160 pm and so 
a disk dominated by large grains. However, there is a high cor- 
relation between a and the pair (J3, Ao) so that these parameters 
cannot be properly constrained in reality. In fact, in the discus- 
sion below, we argue that small grains should be abundant in 
the disk by comparing timescales of collision and removal pro- 
cesses. The best fit values of the other parameters are in Table 
[3] and their uncertainties were determined as described in § 15.21 
The total cross-sectional area of the dust A = 2.3 AU 2 can be 
converted to the dust mass ma -2.2x 10~ 3 y/D max /10 cm in M ffi 
for a collisional cascade, using eq(4) with p = 1.2 g/cm 3 for 
icy grains and D min = 1 pm. The maximum diameter D max is 
unconstrained although objects larger than 10cm contribute neg- 
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Fig. 4. Maps of residuals for the best fit models of the disk to the PACS images at 70, 100, and 160 /urn from left to right. The post-fit 
residuals are between ±3.5<t , coded as black, blue, orange, yellow and white over this range (o"o is the same noise rms as used for 
contours in FigQ]but the color scale is not the same as in FigQ}. In zooming the electronic version, the contours apparent in these 
maps are -3, -2,-1, 1, 2, 3 <x , (dashed contours are negative levels). The two background sources are masked in the 100 and 160 //m 
images as the dotted squares show. 



ligibly to the emission at the wavelengths considered in this pa- 
per. Nonetheless, the size distribution probably extends beyond 
10 cm as discussed in § 18.21 The inclination could be anywhere 
within a relatively broad range (30° < i < 70°) that matches the 
purely geometrical determination based of the ratio of the major 
and minor axes of the Gaussians fit to the central emission in 
§ 14.21 The inner radius is 25 ± 12 AU potentially providing an 
indication of the scale of the planetary system around GJ 58 1 . In 
a similar way to the fits of the individual images, we cannot dis- 
tinguish between a relatively narrow ring and a disk extending 
beyond 100 AU with this combined fit. The best fit value of fj is 
3.5^j g, making the dust temperature between 50 and 30 K over 
the extent of the disk, despite the low luminosity of GJ 58 1 . This 
factor is partially correlated with A as shown in Fig [5] but it is 
clearly inconsistent with unity as we have established by forcing 
fr = 10 in the model and found a reduced xl as high as 1 -92 
for the best fit with this constraint. An emission model including 
a grain size distribution instead of a fixed fj as in our current 
model would be more realistic, but this can be best implemented 
only if the SE P is finely sampled sp ectroscopically from mid-IR 
to submm Ce.g. lLebreton et al.ll2012h . 



5.4. Model with a Gaussian profile for the grain surface 
density 

Finally, instead of a power law for the radial distribution of the 
grain surface density, we tested a Gaussian profile S g exp(-0.5 x 
((r - r g )/w g ) 2 ) peaking at radius r g and having FWHM of 
w g x 2 a/2 ln2. The sky grid for the model, the calculation of 
xl, and the ranges of model parameters tested for A, fr, i, and Q 
were as in § 15.21 Values of r g ranged from 2 X w g to 150AU, and 
w g ranged from 3. 1 AU to r g /2 (the Gaussian profile is truncated 
to 2w g toward the star). We modeled the three PACS images indi- 
vidually and in combination. The resulting best fits have reduced 
xl of 1.03, 0.94, and 1.05 at 70, 100, and 160 fim, and of 1.04 
for the combined images with j3 = (no grey body effect). The 
residual maps are featureless. These best fits are statistically in- 
distinguishable from those with the power law presented in § 15.21 
and !5.3l The resulting parameters are : r g = 52+ 15AU, FWHM= 
38 + 15AU, A = 2.5 ± 1.2 AU 2 , f T = 3.0 + 0.5, i < 60°, and 
a = 120 + 20°. 



In this model, the inner part of the system is populated with 
dust making the inner radius determined with the power law sur- 
face density in § !5.3l less definitive. 



white ; X min 
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Fig. 5. Map of the reduced xl showing the correlation between 
the temperature factor fr and the total cross-sectional area A of 
the dust. The minimum^ 2 is the white region in this map. 



6. The SED and IRS spectrum of GJ 581 

We present the SEDs of the star GJ 581 and modeled dust emis- 
sion that we used to determine the fractional dust luminosity 
Ldust/L* ~ 10~ 4 of the disk. The archival IRS spectrum shows 
a marginally significant excess above the photospheric level that 
provides additional constraints on the dust emission. However, 
we show that a single cold disk model and a two component 
disk model cannot be distinguished to explain this 2-cr excess. 



6. 1 . The SED and the fractional dust luminosity of the disk 

The photometry data collected for the SEDs are summarized 
in Table [4] The flux densities have been color corrected when 
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A Cum) A (/tm) 

Fig. 6. SED of the cold disk model. The Left-hand figure shows the best fit to the three PACS images only (xl = 1 -03, Table|3]l. The 
Right-hand figure shows the best fit to both the three PACS images and the IRS spectrum (ix 2 PACS +^ R5 )/2 = 1.05). The modeled 
cold dust emission is the blue curve, the Next Gen stellar atmosphere spectrum is the grey curve, and their sum is the green curve. 
The IRS spectrum is in red. The upper inset zooms on the IRS wavelengths and displays spectra as S v x (A/24) 2 on a linear scale for 
clarity. The lower inset zooms on the three PACS bands. In the left-hand figure, the best fit model satisfactorily fits the PACS data 
as shown by the lower inset but misses the IRS data as shown by the upper inset. In the right-hand figure, the best fit model partially 
misses the PACS data but satisfactorily fits the IRS data. 



required Q The SED of GJ 581 is based on th e NextGen stel- 
lar atmospheric model dHauschildt et al.l ll 9991) . with the value 
log(g)=5.0 and the effective temperature 3500 K, fit to the 
Johnson UBV and Cousins RJ photometry, the JHK S photom- 
etry from 2MASS, and the recent photometry from AKARI and 
WISE. Note that the flux densities of the photosphere used for 
our modeling in § [5] were predicted from this fit (5.8, 2.8 and 
1.1 mJy at 70, 100 and 160 /mi, respectively). In Fig [6] (left-hand 
panel), we show this SED for the star and the SED for the dust 
emission from our modeling. The fractional dust luminosity was 
determined by integrating the SED of the dust emission and is 
Ldust/L* = 8.9 x 10~ 5 . This value is consistent with the fractional 
dust luminosity Q a i, s . A / An.r 2 = 9.9 X 10" 5 determined from 
the cross-sectional area of the grains A - 2.3 AU 2 from our fit 
in Table [3] using the mean disk radius r = (25 + 60) /2 = 43 AU, 
and assuming the absorption efficiency Q„b s = 1 for the grains 
larger than 1 /mi. The agreement between these two indepen- 
dent determinations of the fractional dust luminosity provides a 
self-consistency check of our modeling. This fractional dust lu- 
minosity is higher than that of the Kuiper Belt by several orders 
of magnitude. 

6.2. IRS spectrum 

6.2.1. Synthetic photometry 

The Spitzer IRS spectrum is superimposed on the star's SED 
in Fig [6] As is standard with IRS spectra, the short wavelength 
module SL (7.6 - 14.2 //m) has to be adjusted to the predicted 
photosphere, and IRS flux densities were scaled up by the factor 
1.066 for GJ 581. In Fig|6]and insets, a small excess is apparent 
above the photospheric level at the longest wavelengths of the 
spectrum (module LL1 : 20.4 - 34.9 /mi). 

We have carried out synthetic photometry with a rectan- 
gular bandpass between 30 and 34 //m and between 15 and 
17 /mi which gives the widest wavelength range while still inside 

7 http://herschel.esac.esa.int/twiki/pub/Public/PacsCalibrationWeb/ 
cc_report_vl.pdf 



of the Long -Low IRS module as prescribed in Carp enter et al.l 
(2008, 2009). We computed the synthetic flux densities Sn.^m = 
32.3 ± 1.9 mJy (IRS) and 28.4 mJy (Next Gen) yielding the 2<x 
excess 3.9 ± 1.9 mJy, and S i 5 .9 6/ ™ = 110.7 ± 0.85 mJy (IRS) 
and 109.2 mJy (Next Gen) yielding the lower significance ex- 
cess 1 .5 + 0.85 mJy. We computed these synthetic flux densities 
as the weighted mean of the data points in these bands, and using 
the same weights for the corresponding Next Gen synthetic flux 
densities. The IRS flux density uncertainty includes an absolute 
calibration error of 6%. Photospheric flux densities predicted for 
late type stars (K and M) by the Kurucz or Next Gen models 
have been shown to be overestimated in the mid -IR by as much 
as 3-5% dGautier et al.l2007HLawler et al.l2009b . Hence, the sig- 
nificance of the marginal excess at 3 1 .6 /mi is likely higher in re- 
ality. If real, this excess for the mature M-star GJ 581 is notable 
because, even among A-type and solar-type stars, 24 /mi ex- 
cesse s are less frequen t than 70 /mi excesses and decrease wit h 
age dRiekeet al.1 120051 iTrilling et al.(l2008b lLohne et al.1 12008). 
In the next two sections, we investigate the implications for the 
system around GJ 581 if this excess is real. 

6.2.2. Modeling the IRS and PACS data with the cold disk 
model 

First, we fit the single cold disk model of § [5] simultaneously 
to the three PACS images and the IRS spectrum, minimizing 
Xm = (xIacs + X 2 ir S )I 2 where *p AC s and *«s are the reduced 
xl for the PACS and IRS data, respectively. With this definition, 
both data sets have the same weight in the fit. The best fit model 
thus obtained is characterized by xfoi - 1-18, resulting from 
X 2 IRS = 1.25 andx 2 PACS = 1.11, and its SED is shown in Fig [6] 
(right-hand panel). The main parameter changes are fj =5.5 
and A = 0.8 AU 2 , instead of f T = 3.5 and A = 2.3 AU 2 in 
Table [3] This value ofx 2 or is higher than^ 2 = 1 .03 of the best fit 
in this Table and is high for the number of degrees of freedom 
of 1186 in ^-statistics (probability = 1% of pure noise). It is 
instructive to compare the SEDs of these two fits in Figure [6] ; 
the simultaneous fit to the PACS and IRS data in the right-hand 
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A (,um) 

Fig. 7. SED of the best fits of the cold disk and warm belt mod- 
els. The green curve corresponds to the warm and cold dust 
emissions added to the Next Gen stellar atmosphere spectrum 
{grey). IRS spectrum is in red. We also show separately the 
best fit of the warm belt emission (r w = 0.2 AU and m w - 
2.8x1 (T 6 M(j ) (yellow) and the best fit of the cold disk emission 
(parameters of the combined fit are in Table 0) (blue). Insets and 
photometric data points are the same as in the legend of Fig [6] 



panel does appear to be skewed to some degree. The assump- 
tion in our current model that the temperature does not depend 
on grain size and wavelength is a limitation. A size distribution 
would broaden the SED, and it may improve the ability to fit a 
single disk model to the flux densities of the IRS spectrum and 
the PACS bands simultaneously. 

6.2.3. Modeling the IRS and PACS data with a two 
component model 

We explore another possibility, a two component model in which 
a belt of warm dust is added to our cold disk model of |5] The 
model of this belt is simply based on blackbody grains (i.e., we 
set fj — 1 for the warm component) located at radius r w and 
having a total cross-sectional area A w . This two parameter model 
is fit to the IRS spectrum alone by minimizing x^ RS . We found 
r w = 0.2 AU (T dusl = 191 K) and A w = 7 x 10 5 AU 2 , giving 
a corresponding dust mass of m w = 2.8 x 10~ 6 M(j , assum- 
ing the standard grain size distribution (oc D~ 35 ) between 1 fim 
and 1 mm-sized particles and p = 3 g/cm 3 (noting the depen- 
dence of this estimate on the unknown maximum size given in 
eq. 4). However, acceptable fits could also be found for r w be- 
tween 0.05 AU (T dust = 382 K) and 0.4 AU (T dust = 135 K), en- 
compassing the orbital radii of the planets GJ 58 lc and GJ 58 Id. 
The IRS data alone cannot constrain fj but if this parameter were 
larger than unity, the dust would be fj times further out than the 
r w quoted above for the corresponding T dust . The SED of the 
two component model is shown in Fig [7] where we had to de- 
crease the cold dust cross-sectional area A by 6% from its value 
of Table[3]to account for the warm dust contributions to the flux 
densities at 70, 100 and 160 //m. The fractional dust luminosity 
is L dust /L* = 5.7 x 10" 5 for this warm belt shown in yellow in 
Fig [7] Such a belt is comparable to the warm disk around the 
K0 star HD69830 dLisse et al]l2007l) . However, the proximity of 
the warm dust to the kno wn planets suggests that it could be dy- 
namically unstable (e.g. lMoro-Martm et al.1 12007*). Nevertheless, 



definitive proofs are still missing to establish the reality of this 
warm belt in the GJ 581 system. 



Table 4. Photometry of GJ 581 



Wavelength 


s ,. 


References 


(yum) 


(mjy) 




0.36 


8.3 ±2 


Hipparcos (Koen et al. 2010) 


0.44 


61 ± 10 


(") 


0.55 


222.7 ± 3 


(") 


0.66 


523.1 ± 13 


(") 


0.81 


1490 ± 14 


(") 


1.23 


3317 ± 82 


2MASS fCutri et al. 2003) 


1.66 


3939 ± 120 


(") 


2.16 


3051 ± 65 


(") 


9.0 


322 ± 18 


AKARI (Ishihara et al. 2010) 


11.6 


213 ± 19 


WISE (Wrieht et al. 2010) 


22.1 


61.2 ±6 


(") 


70.0 


18.9 ± 1.4 


PACS this work 


71.42 


20.0 ±5.3 


MIPS (*) 


100.0 


21.5 + 1.5 


PACS this work 


160.0 


22.2 ± 5.0 


(") 


250.0 


< 24° 


SPIRE this work 


350.0 


< 26 (3<r) 


(") 


500.0 


< 27 (3a-) 


(") 


1200.0 


< 2.1 (3cr) 


MAMBO (Lestrade et al. 2009) 



Color correction factors : 1.125 for AKARI, 0.956 for WISE at 
1 1.6 jum, 0.987 for WISE at 22.1 /im, and 0.992, 0.980, 0.995 for PACS 
(T=40 K) at 70.0, 100.0, 160.0 ^m, respectively. 
(*) : not used in the fit. 
a : see § 1441 



7. Brightness limit on scattered light around GJ 581 

Our HST/NICMOS Fl 10W image of GJ 581 after PSF subtrac- 
tion (Fig. [8]) is sensitive to a region from 4" radius (30 AU) 
to approximately 10" radius (62 AU) along PA=120 degrees. 
We estimate the 3cr sensitivity to nebulosity in this region as 
Efiiow = 1 8.7 mag arcsec' 2 . We used the radiative transfer code 
MCFOST (iPinteet all [20061) to produce synthetic F110W de- 
bris disk images for pure astronomical silicate and pure wa- 
ter ice models that match the SED based on the geometry de- 
rived in Table [3] using the standard F110W NICMOS through- 
put. The maximum surface brightnesses in the 4 - 10" range 
at the forward scattering peak for pure water ice grains and 
for pure astronomical silicates are ~ 19.4 mag arcsec~ 2 and 
~ 21.2 mag arcsec~ 2 , respectively. Hence, our dust models 
are consistent with the non-detection of scattered light around 
GJ 581, but show that the disk may be detectable in deeper ob- 
servations. 



8. Discussion 

We have spatially resolved a disk around the mature M-star 
GJ 581 hosting four planets. This cold disk is reminiscent of 
the Kuiper Belt in the Solar system but it surrounds a low mass 
star (0.3 M ) and has a much higher fractional dust luminosity 
Ldust/L* of ~ 10" 4 . It shows that debris disks can survive around 
M-stars beyond the first tens of Myr after the protoplanetary disk 
disperses, and they can be detectable although they have been 
elusive in searches so far. 
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Fig. 8. HST/NICMOS F110W image of GJ 581 (North is up, 
East is left). The star was not placed behind the occulting spot 
available on the NIC2 camera. We subtracted the PSF using 
observations of GJ 250 B made ear l ier in the same scientific 
program, as described in iKrist et al.l (Q~998). The circular black 
digital mask has 4" radius (30 AU) and blocks the central re- 
gion where PSF subtraction artifacts are significant. Along the 
position angle of the disk (PA=120 degrees) the field of view 
is limited to approximately 10" radius (62 AU). We estimate 
the 3<x sensitivity to nebulosity in the 4"-10" radius region as 
Sfiiow = 18.7 mag arcsec~ 2 . Lack of detectable scattered light 
at this level is consistent with the dust model derived from the 
far-IR PACS images. 



8. 1. Dust temperature in the cold disk 

The factor fj is significantly larger than unity in our analysis and 
indicates that the dust temperature ranges from ~ 50 to ~ 30 K 
from the inner to the outer radius of our modeled disk. This is 
about three times the black body equilibrium temperature for 
the dust around this low luminosity M3-type star. Values of fj 
larger than unity have als o been found for th e debris disks around 
the G -type star 61 Vir dWvatt et alj 1201 2|) and several A-type 
stars (Boo th etal.ll20H . This is akin to disks resolved in scat- 
tered light which tend to be m ore extended than their sizes es - 
timated from blackbody SED (Rodrig uez & Zuckermanl l2012). 
Values fj > 1 are interpreted as evidence for dust grains of small 
sizes and/or optical propert i es different from b l ackbody spheres 
dBackman & Parescd [19931: iLisse et all 120071: iBonsor & Wvattl 
2010). When the SED of a debris disk is finely sampled spectro- 
scopically from mid-IR to submm, a parameter search for com- 
position, structure, si ze distribution of the grains can be con- 
ducted usefully (e.g. Lebre ton et alj I20T2I) . Such a parameter 
search would be degenerate for GJ 581 because of its limited 
photometry, and so these effects have been reduced to the single 
parameter f T in our model. 

8.2. Collision, Poynting-Robertson and stellar wind time 
scales for the GJ 581 system 

In addition to gravitational forces, dust dynamics is controlled 
by radial forces (radiation and stellar wind pressures) and by 
tangential forces (Poynting-Robertson and stellar wind drags). 
For large dust grains, these perturbing forces act on much longer 
timescales than collisions, and such grains simply orbit the star 
until they are broken into smaller fragments in collisions with 



other grains. This results in a collisional cascade with a size 
distribution with a characteristic slope n(D) oc D~ 7 ^ 2 (assum- 
ing dust grain strength is independent of size). For small dust 
grains, perturbing forces truncate (or at least significantly de- 
plete) the size distribution at scales where one of the perturb- 
ing force timesca les is shorter than the collisional lifetime (e.g. 
IWvatt et alJl20TTh . To ascertain the process dominating the dust 
removal requires a comparison of the relevant timescales. 

Fig. [9] shows the ratio (3 of the radiation pressure to stellar 
gravity experienced by icy dust grains of different sizes. This 
peaks for sizes comparable to the wavelength where the stellar 
spectrum peaks and is also proportional to L,/M,,, where L , and 
M» are the luminosity and mass of the star (lGustafsonll 994). The 
low luminosity of GJ 581 means that f3 < 0.5 for all icy grains 
regardless of size, and the same is true for other compositions. 
Since dust with (3 < 0.5 that is created in collisions is always 
placed on a bound orbit, this means that radiation pressure is 
not a m echanism that can be invoked to expel the dust from the 
system dWvatt et al.lll999t) . 

Fig. [9] also sh ows the ratio fi sw of the stellar wind pressure 
to stellar gravity dGustaf son! 1 1 9941) . This depends on the stel- 
lar mass loss rate and stellar wind speed that are poorly under- 
stood and hard to measure for M-stars. Here we estimate the 
mass loss rate from the non-detection of X-rays from GJ 58 1 
by ROSAT implying log L x < 26.44 erg/s dSchmittet al.lll995l) . 
The correlation b etween X-ray surf ace flux and mass loss rate of 
GKM-type stars dWood et alJl2005h then yields the upper limit 
of 2 M , where the solar mass loss rate M G = 2 x 10~ 14 M yr _1 . 
We also consider the stellar w ind speed to be ~400 km/s as ap- 
propriate for GKM-type stars (Wood 2004). Finally we assumed 
100% efficiency of momentum coupling b etwee n dust and the 
stellar wind (and used eq. 12 of lPlavchan et al.1 d2005l) ). With 
these assumptions we found that f3 sw can only be > 0.5 for 
dust smaller than a few nm, meaning that stellar wind pressure 
could only truncate the collisional cascade below the nm-scale; 
furthermore, stellar wind pressure would be ineffe ctive if small 
grain s couple inefficiently to the stellar wind (e.g. Minato et al. 
12009) . 

A comparison of timescales first requires an estimation of 
the collisional lifetime of dust grains of different sizes. Here we 
use eq. 4 and the parameters for the disk found from the mod- 
eling of the combined images presented in Table [3] to derive the 
total mass M tot = 2.2 x 10~ 3 y/D c /lQ cm in M®, assuming the 
standard size distribution between D m! „ = 1 /im and the diameter 
D c of the largest objects and the density p = 1.2 g/cm 3 f or icy 
grains . The collisional lifetime is estimated using eq. 16 o flWvatj 
(2008) with the additional assumptions that dust orbital eccen- 
tricities are 0.05 and that their strength is 10 3 J/kg (independent 
of size so as to be consistent with the assumptions about the size 
distribution). The resulting collisional lifetime is 0.22 V75 Myr, 
where D is in /im. Since we expect the cascade to extend up to 
sizes for which their collisional lifetime is equal to the age of the 
star (~ 5000 Myr), then as long as our assumptions apply up to 
large sizes we can get a rough estimate of the total mass of the 
disk as 0.16M e in objects up to D c = 0.5 km in diameter. 

Fig. [TO] shows the timescale for dust to migrate from the in- 
ner edge of the disk at 25 AU to the star due to P-R drag. The 
dependence of this timescale on particle size results from a scal- 
ing oc l/p which means that this has a minimum value of 60 Myr. 
Since this timescale is longer than the collisional lifetime at all 
sizes, P-R drag is not a significant loss process from the disc. 
Fig. [TO] also shows the corresponding timescale for migration 
due to stellar wind drag. This also includes a scaling oc l//3 sw , 
and the efficient momentum coupling assumed here means that 
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Fig. 9. The ratios of radiation pressure (solid line) and stellar 
wind pressure (dashed line) to stellar gravity as a function of 
particle diameter for icy grains around GJ 58 1 . Particles with /3 > 
0.5 are put on hyperbolic orbits as soon as they are created and 
so removed from the system on orbital timescales, thus setting 
the blow-out limit. 

this timescale decreases indefinitely to smaller sizes oc D. As a 
result, stellar wind drag timescales become shorter than colli- 
sional timescales at a size of around 3 nm. 

Thus, if all of the assumptions hold, we would expect the 
collisional cascade to extend down to 3 nm, while smaller dust 
is removed by stellar wind drag. However, it should be noted that 
there are significant uncertainties, both on the magnitude of stel- 
lar wind drag and its efficiency of coupling to small grains, and 
on the geometry of the dust disk which impacts the collisional 
lifetimes. As such this plot should be considered as representa- 
tive of the kind of arguments that need to be considered when as- 
sessing the fate of material in the debris disk of GJ 581. Further 
study of this issue is left to later papers, but here we note that the 
existence (or not) of grains smaller than 1 fim is not important 
for the observable properties of the disk discussed in this paper, 
since such grains are inefficient emitters in the far-IR. 

Another scenario that we have not considered in detail is 
that the dust is all in large mono-sized grains, in a configura- 
tion meaning that the dust collides at low enough velocities that 
particles bounce off each other rather than destroy each other 
dHeng & Tremainell2010l) . Two constraints on such models are 
that the SED should look like a black body (since all the dust is 
large), and the fractional luminosity should not be large enough 
that collisions must necessarily occur at high velocity. Here the 
fractional luminosity only constrains the collision velocity at the 
inner edge to be > 0.3m/s which is not sufficient to require a col- 
lisional cascade. However, although there is no evidence from 
our limited photometry of GJ 58 1 that the spectrum departs from 
black body shape, the resolved location of the dust shows that it 
is significantly hotter than black body, consistent with the pres- 
ence of small grains and so incompatible with this model. 

8.3. Planets and disk relationship for the GJ 581 system 

First, we note that our determination of the inclination of the 
disk relative to the plane of the sky is 30° < i < 70° (face on 
disk is i = 0°). This is mostly constrained by the 160 //m image 
and is fairly insensitive to the masks used for the background 



Fig. 10. Dust removal timescales as a function of particle size, 
due to collisions (solid line), Poynting-Robertson drag (dashed 
line), stellar wind drag (dash-dot line), and stellar wind pressure 
(dotted line). 



sources. If the disk mid-plane and the orbits of the planets are 
coplanar, this range of inclination makes the masses of the plan- 
ets of GJ 58 1 no more than ~ 1 .6 times their measured minimum 
masses by radial velocity and, interestingly, ensures the long- 
term stabi lity of the orbits in this system as shown in dynamical 
studies bv lBeust et ail d2008l) and lMavor et all d2009h . 

In our DEBRIS sample of 89 M-stars, there are only three M- 
stars with known planets overall (GJ 876, GJ 832 and GJ 581). 
GJ 581 hosts low mass planets and now has a detected disk, 
while GJ 876 and GJ 832 host Jupiter mass planets and have no 
detected disk brighter than the fractional dust luminosity 10~ 5 in 
our survey as we shall present in a future study (Matthews et al. 
in prep). Hence, using these three stars as a sample, the outcome 
is one disk for one low mass planet system (1/1) and no disk for 
two high mass planet systems (0/2). Although this is small num- 
ber statistics, we note that it is suggestive that the correlation 
between l ow-mass planets an d debris disks recently found for 
G-stars by IWvatt et alJ ([2012) also applies to M-stars. It is also 
intriguing that the only debris disk confidently detected in our 
current analysis surrounds the one star in the sa mple that hosts 
low-mass pl anets. We note that simulations by Ray mond et al.l 
d20ni2012h suggest that a correlation might exist between low 
mass planets and debris as a result of planet formation processes. 
Note that the star AU Mic does not fall in the DEBRIS sample 
of the nearest M-stars, and so is not included in the statistics 
above; this young star (12 Myr) has a bright disk, but no known 
planets, although radial velocity measurements toward AU Mic 
are insensitive to planets with masses lower than a few Jupiters 
even for short orbital periods bec ause of its high ch romospheric 
activity (see GJ 803 in Fig 19 of lBonfils et alfeoilh . 

Current programs show that a large fraction of M-stars are 
orbited by low-mass p lanets. The radial ve locity survey of 102 
M-stars conducted by iBonfils et alJ d201 ll) yields the high oc- 
curence of 35^|j% for low mass planets (2-10 M ffi ) around 
M-stars, unlike the low occurence of giant planets of ~ 2%, 
for orbital periods under 100 days. Transit observations in the 
Kepler field show that small candidate planets (2 - 4R @ ) with 
P < 50 days are found around 25 + 10 % of the M-stars 
(T e ff = 3600 - 4100 K), seven times more frequently than 
around the hottest stars (6600-7100 K). There is no such a de- 
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pendence for larger planets (4 - 32R e ) with P < 50 days, found 
uniformly arou nd 2 + 1 % of the star s across all spectral types in 
the Kepler field dHoward et al.l2012b . Hence, if there is a correla- 
tion between the debris disks and low-mass planets of M-stars at 
a level similar to that found for G stars (4/6 nearby G- stars with 
detecta ble low-mass planets have detectable disks, I Wvatt et ail 
(1201 2l) . the high fraction of M-stars with low-mass planets would 
explain the detection of the disk around GJ 581 but would imply 
also that more disks were expected to be detected in our DEBRIS 
M-star sample. We defer this discussion to a paper that describes 
those observations in more detail, since not all observations are 
sensitive to disks at the same level. However, GJ 58 1 is at the me- 
dian distance of our DEBRIS M-star sample that was observed 
to uniform depth, so it could simply be the brightest because of 
some intrinsic properties. The explanation could also be related 
to its multiple planetary system. 

Secular perturbation theory has been a pplied to planetesi- 
mals in debris disks perturbe d by planets in I Wvatt et al.l (l999) 
and IMustill & Wvattl (120091) . The outermost planet GJ 581d 
(5.4 M e , a p /=0.22 AU, and e p /=0.25, the highest eccentricity in 
the system) cannot stir the disk at a — 25 AU because the time 
scale for orbital crossin g of planetesimals is m uch longer than 
the stellar age (eq. 15 in Mus till & Wvattll2009l) . However, a hy- 
pothetical outerplanet, for example a Neptune mass planet (17 e ) 
at 5 AU with a moderate orbital eccentricity of 0.2, can stir the 
disk at a = 25 AU in much less than the age of the system, and 
trigger destructive col lisions of 0.5 km-sized bodies (eq. 27 in 
Mustill & Wvattl 12009b to feed a collisional cascade. The most 
recent detection limit on m sini from radial velocity measure- 
ments of GJ 58 1 over 3. 3 years indicates tha t such a planet would 
not have been detected (Bonf ils et al.ll201 ll Fig 13), and there is 
a large region of parameter space of m sini vs a over which a 
planet could both stir the disk and have eluded detection in ra- 
dial velocity measurements. 

Alternatively, mild collisions between planetesimals in 
a weakly excited disk could eventually form a Pluto-sized 
body, which in turn stirs the disk so that it produces dust 
dKenvon & Bromlevll2008l) . This self-stirring scenario is plau- 
sible since the time scale required for the formation of a Pluto- 
sized body at 50AU is comparable to the age of GJ 581, even 
when the surface density of solids is ten times smaller than 
the minumim-mass solar nebula around an M3-type star (eq. 41 
iKenvon & Bromlevll2008L with x m = 0.1 ). 

9. Conclusion 

We have spatially resolved a debris disk around the M-star 
GJ 581 with Herschel PACS images at 70, 100 and 160 fxm 
and modeled these observations. This is the second spatially re- 
solved debris disk found around an M-star after AU Mic, but, 
in contrast, GJ 581 is much older and is X-ray quiet. Our best 
fit model is a disk, extending radially from 25+12 AU to more 
than 60 AU. Such a cold disk is reminiscent of the Kuiper Belt 
but it surrounds a low mass star (0.3 M Q ) and its fractional dust 
luminosity L^ ust IL„ of ~ 10~ 4 is much higher. Also, in our best 
fit model, the dust temperature is found to be significantly higher 
than the blackbody equilibrium temperature indicating that small 
grains are abundant. Finally, the inclination limits of the disk 
make the masses of the planets small enough to ensure the long- 
ter m stability of the s yste m according t o dyn amical simulations 
bv lBeust et all (120081) and lMavor etall (120091) . 

This disk complements our view of this remarkable system 
known to host at least four low mass, close-in planets. These 
planets cannot perturb sufficiently the modeled cold disk to trig- 



ger destructive collisions between planetesimals over the age of 
the star, but a hypothetical outer planet, for example a Neptune 
mass planet with an orbital radius of 5 AU and a moderate ec- 
centricity, could replenish the system with dust. Alternatively, 
the self-stirring mechanism could operate for this old star caus- 
ing sufficient dynamical excitation to produce the observed dust. 

It is intriguing that, in our current analysis of the DEBRIS 
sample of 89 M-stars, the only debris disk confidently detected 
around a mature M-star also happens to be around the only star 
known to have low mass planets. This could mean that the corre- 
lation between low-mass p l anets and debris disks recently found 
for G-stars by Wva tt et all d2012l) also applies to M-stars. Then, 
the high fraction (~ 25%) of M-stars known to host low mass 
planets in the radial velocity and Kepler observations should 
make debris disks relatively common around them. If these disks 
have not been detected yet, it may be because searches have sim- 
ply not been deep enough, or because the disk around GJ 581 is 
the brightest owing to some intrinsic properties ; for example 
hosting a multiple planetary system. 

Future studies and complementary observations of GJ 581 at 
higher angular resolution will enhance further our knowledge of 
this remarkable system around an M-star. 
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